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STREAMFUNCTION FORMULATION OF THE STATIONARY 
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Abstract. In this paper we propose a two- level finite element discretization of the nonlinear 
stationary quasi-geostrophic equations, which model the wind driven large scale ocean circulation. 
Optimal error estimates for the two-level finite element discretization were derived. Numerical experi- 
ments for the two-level algorithm with the Argyris finite element were also carried out. The numerical 
results verified the theoretical error estimates and showed that, for the appropriate scaling between 
the coarse and fine mesh sizes, the two-level algorithm significantly decreases the computational time 
of the standard one-level algorithm. 
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1. Introduction. Two-level algorithms are computationally efficient approaches 
for finite element (FE) discretizations of nonlinear partial differential equations (3) [4] 
ITT1 Uni [25J ■ A two-level FE discretization aims to solve a particular nonlinear elliptic 
equation by first solving the nonlinear system on a coarse mesh and then using the 
coarse mesh solution to solve the linearized system on a fine mesh. The apeal of such a 
method is clear; one need only solve the nonlinear equations on a coarse mesh and then 
use this solution to solve on a fine mesh, thereby reducing computational time without 
sacrificing solution accuracy. The development of the two-level FE discretization was 
originally performed by Xu in [55]. Later algorithms were developed for the Navier- 
Stokes equations (NSE) by Layton Q35] (see also [H] IT2 ] 1251 1201 1251 121] ) and for the 
Boussinesq equations by Lenferink |22) . 

As computational power increases complex models are becoming more and more 
popular for the numerical simulation of oceanic and atmospheric flows. Computa- 
tional efficiency, however, remains an important consideration for geophysical flows 
in which long time integration is needed. Thus, simplified mathematical models are 
central to the numerical simulation of such flows. For example, the quasi-geostrophic 
equations (QGE), a standard mathematical model for wind driven large scale oceanic 
and atmospheric flows [23] [27] , are often used in climate modeling 10] . 

Most FE discretizations of the QGE are for the streamfunction-vorticity formu- 
lation. The reason is that the streamfunction-vorticity formulation allows the use of 
low order (C°) FEs, although one needs to discretize two flow variables, the poten- 
tial vorticity q and the streamfunction ip. We note that the streamfunction-vorticity 
formulation is often used in the numerical discretization of the 2D NSE, to which 
the QGE are similar. Alternatively, one can, instead, use the pure streamfunction 
formulation of the QGE. The advantage lies in an equation that contains only one 
flow variable, the streamfunction, ip, at the price of having to deal with a fourth-order 
partial differential equation. Thus, its numerical discretization with conforming FEs 
requires the use of high-order (C 1 ) FEs, e.g. the Argyris element pQ E] . 

The streamfunction formulation of the QGE still suffers from having to solve 
a large nonlinear system of equations. This is usually done by using a nonlinear 
solver, such as Newton's method. These nonlinear solvers typically require solving 
large linear systems multiple times to obtain the solution to the nonlinear system. 
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Solving these large linear systems multiple times can be time consuming. Thus, a 
two-level algorithm greatly reduces computational time over the standard nonlinear 
solver, since we need only solve the nonlinear system on a coarse mesh and then use 
that solution to solve a linear system on a fine mesh. 

In this paper, we propose a two-level algorithm for the FE discretization of the 
streamfunction form of the stationary QGE (SQGE). This conforming FE discretiza- 
tion is based on the Argyris element. Additionally, we present a rigorous error analysis 
for the two-level FE discretization. The theoretical error bounds as well as the in- 
creased computational efficiency are illustrated numerically for a test problem. 



The rest of the paper is organized as follows: In Section 2 we present the SQGE. 
Then in |Scction 3| we present the weak formulation of the SQGE, including notation 
and functional spaces. Next, |Section 4 contains the presentation of the one-level FE 
discretization of the SQGE. In Section 5 we discuss both the two-level algorithm and 



its application to the SQGE. Next, in Section 6| we provide rigorous error bounds 



for the two-level FE discretization of the SQGE and we discuss the scaling between 



the fine mesh, h, and the coarse mesh, H . Then Section 7 includes numerical results 
which both verify the theoretical error bounds presented in |Section"6| and illustrate 
the computational efficiency of the two-level algorithm over the standard one-level 
method. Finally, in |Scction 8| we present our conclusions. 

2. Streamfunction Formulation. The pure streamfunction formulation of the 
SQGE is 



where 



J(u, v) = 
U 



Ro 



ox 



du dv du dv 
dx dy dy dx 1 
UL 



in f2, 



(2.1) 



(2.2) 



Re" 1 = — 



are the Jacobian, Rosby number, Reynolds number, respectively, and (3, A, U, and L 
are the cooeficient in the beta plane approximation, the eddy viscosity, the charac- 
teristic velocity scale, and the characteristic length scale, respectively (see [T4l 125] ). 
To completely specify (2.1 ), we need to impose boundary conditions (see }9l 1271124] 



for a careful discussion of this issue). In this report, we consider 



dip 
dn 



= on il, (2.3) 
These are also the boundary 



where n represents the outward unit normal to Q. 
conditions used in Q71 D3J [T3] for the 2D NSE - 

3. Weak Formulation. Now we can derive the weak formulation of the SQGE 
in streamfunction formulation (2.1). To this end, we first introduce the appropriate 
functional setting. Let 



dip 
dn 



= on dn 



Multiplying (2.1 1 by a test function and using the divergence theorem, we get, 



in the standard way (see [17] ). the weak formulation of the SQGE in streamfunction 
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formulation: 



where 



Find ip £ X such that 
a(ip, X ) + b(ip] V>, x) + c(V>, X) = *(x), V X £ X, 



a(ip,x) = ReT 1 / AtpAxdx 



(3.1) 



&(C;V>,x) = / A C (V^Xx - V'xXy) 
in 

c(V') X) = -Ro^ 1 / ipxxdx, 



(3.2) 



Jo 



Lemma 1. Given ip, £, ip £ Hq{Q) and F £ H~ 2 (fl), the linear form £, the bilinear 
forms a and c, and the trilinear form b are continuous: there exist T± > and T 2 > 
such that 



a(iP, X ) < Re-^lxU 

&(C;^x)<ri|ci2MaMa 
c{i),x)<Ro- 1 r 2 \i>\ 2 \ x \2 

£(x) < Ro^\\F\\^\x\2- 



(3.3) 
(3.4) 
(3.5) 
(3.6) 



For a proof see [7J . 

For small enough data, one can use the same type of arguments as in [T51 Qj5] 
to prove that the SQGE in streamfunction formulation (2.1) is well-posed. In what 



follows we will always assume that the small data condition involving Re, Ro, and F, 



is satisfied and, thus, that there exists a unique solution ip to (2.1). 

The following stability estimate was proven in |14j : 
Lemma 2. The solution ip of (2.1) satisifies the following stability estimate: 



\ip\\ 2 < ReRo^WFW 



(3.7) 



4. Finite Element Formulation. Let T H denote a FE triangulation of ft with 
meshsize (maximum triangle diameter) H. We consider a conforming FE discretiza- 
tion of ( |3"7lj ), i.e., I fl cl = Hoify. 

The FE discretization of the SQGE dO) reads: Find i\) R £ X H such that 



a^ H , X H ) + b^ H ;^ H ,x H )+c^ H ,x H )=e(x H ), V X ff G X H . (4.1) 

Using standard arguments [TIjJ [TC], one can prove that, if the small data condition 
used in proving the well-posedness result for the continuous case holds, then (4.1) 



has a unique solution i/j h (see Theorem 2.1 in [7]). Furthermore, one can prove the 
following stability result for ip H using the same arguments as those used in the proof 
of Lemma [2] for the continuous setting (see [T4"]). 

Lemma 3. The solution ip H of (4.1) satisfies the following stability estimate: 



\^ H \ 2 < ReRo- 1 \\F\\^ 
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(4.2) 



As noted in Section 6.1 in |8] (see also Section 13.2 in |T7|, Section 3.1 in fl8l 



and Theorem 5.2 in [3]), in order to develop a conforming FE for the SQGE (3.1), 
we are faced with the problem of constructing subspaces of Hq(H). Since the stan- 
dard, piecewise polynomial FE spaces are locally regular, this construction amounts 
in practice to finding FE spaces X H that satisfy the inclusion X H C C 1 (57), i.e., find- 
ing C 1 FEs. Several FEs meet this requirement (see, e.g., Section 6.1 in [8], Section 
13.2 in [T7], and Section 5 in [S]): the Argyris triangular element, the Bell triangu- 
lar element, the Hsieh-Clough- Tocher triangular element (a macroelement), and the 
Bogner-Fox-Schmit rectangular element. In this study we use the Argyris element. 

5. Two-Level Algorithm. In this section we propose a two-level FE discretiza- 



tion of the SQGE ( pTlj ). We let X h , X H C X = H$(Q) denote two conforming FE 
meshes with H > h. The two-level algorithm consists of two steps. In the first step, 
the nonlinear system is solved on a coarse mesh, with mesh size H . In the second 
step, the nonlinear system is linearized around the approximation found in the first 
step, and the resulting linear system is solved on the fine mesh, with mesh size h. 
This procedure is as follows: 

Algorithm 1 Two-Level algorithm 

Step 1: Solve the following nonlinear system on a coarse mesh for ip G X H : 

a^ H lX H )+b(4> H ;i> H ,x H ) + c^ H ,X H )=e(x H ), for all X H G X H . (5.1) 
Step 2: Solve the following linear system on a fine mesh for ip h G X h : 

a^'\ X h )+b^ H ^\x h ) + c^\x h )=e(x h ), for all X h € X h . (5.2) 



The well-posedness of the nonlinear system was proven in [7] , see also [14] . The 
following lemma proves the well-posedness of the linear system ( |5.2[ ). 

The following theorem provides an error estimate for the approximation in Step 
1 of the two- level algorithm ( Algorithm 1 ) . 



Theorem 1. Let i/j be the solution of (3.1) and ip H be the solution of (5.1). Fur- 
thermore, assume that the following small data condition is satisfied: 

Re- 2 Ro > ri||F||_ 2 . 

Then the following error estimate holds: 



IV - < C(Re, Ro, r l7 T 2 , F) inf |</> - X H | 2 , 

x H ex H 



(5.3) 



whe 



C(Re,Ro,T 1 ,r 2 ,F) 



r, Ro- 



2Re- 



r 1 i?ei?o- 1 ||F||_ ; 



Re- 



Ti ReRo- 1 ||F||_2 



For a proof see [14] . 

Lemma 4. Given a solution ijj H of (5.1), the solution, ip , to (5.2) exists uniquely. 
Proof. First we introduce the bilinear form B : X h x X h - 



I given by 

5(VA x h ) = a(i> h , X h ) + b(^ H ; 4>\ X h ) + <iP h , X h )- 
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(5.4) 



By [Lemma 1| 



X h ) < {Re' 1 + T^ H \ 2 + Ro- 1 T 2 ) |^| 2 \ X %, W> h , X h € X> 



(5.5) 



The stability estimate, for Tp H , in Lemma ~3| and inequality (5.5 1 imply that B is 
continuous. Additionally, the fact that b(ip H ; ip h , = and c(ifj h ,?p h ) = for all 
?/>' 1 G A' 1 combined with the Poincare-Friedrichs inequality gives 

B(ip h , ip h ) > C \\ip h \\2, y^ h ex h . 

Thus, B is coercive. Therefore, by Lax-Milgram, ip h exists and is unique. □ 
In addition to the existence, uniqueness, and stability of the solution to the con- 



tinuous linear system (Lemma 4) we also have a stability bound for the solution on 



the discrete fine mesh, h. 

Lemma 5. The solution, ip , to (5.2) satisfies the following stability bound: 



U h h<ReRo- l \\F\\_ 2 , (5.6) 

Proof. Setting \ h = ^ m ( |5.2[ ), and noting that c(tp h ,x h ) = — c(x'\ "tp 1 *), which 
implies that c(ip h ,ip h ) — 0, and b(ip H ; ip h , ip h ) — gives 



\1>% = Re^r < ReRo-HFW-a, 



where in the last inequality we used (3.6). Therefore, it follows that 



||^|| 2 <i?ei? - 1 ||^||_ 2 . □ 
6. Error Bounds. The main goal of this section is to develop a rigorous nu- 



merical analysis for the two-level algorithm ( Algorithm 1 ) . The proof for the error 
bounds follows the pattern presented in 

We first introduce an improved bound on the trilinear form 6(C;£,x)- To this 
end, we use the following discrete Sobolev inequality 



|V¥> h || £ ~ <cJ\MK)\\y% y v h ex h . 



(6.1) 



The following lemma will be useful in determining the error bounds for Step 2 
of the two-level algorithm. The first lemma which corresponds to Lemma 5.1 in |llj . 
follows from (6.1 1 and (3.4) and places error bounds on the trilinear form b(tp; x > £)• 
Lemma 6. For any x h € X , the following inequality holds: 



M^x h -X)\<cV\^W\m^\i\x h \2, 



whe 



M£; XA)= / [XyCxy - ZxXyyWy - iiyXyx ~ ZyXxxtyx dx. 



(6.2) 



The following lemma, which corresponds to Lemma 5.6 in |llj . will be useful for 
proving the error bounds for | Algorithm 1[ by allowing one to permute the terms of 
the trilinear term: 



Lemma 7. For ip, £, x G we have 

&(V>; e, x) = &o(£; x, VO - &o(x; e, V0- (6.3) 

The following theorem gives the error bound after Step 2 of the two-level algorithm 
(Algorithm 1 ) and is the main result of this paper. The proof of this theorem is similar 
to the proof of Theorem 5.2 in [TT] . 

Theorem 2. Let ip be the solution to (3.1| and ip h the solution to (5.2). Then ip h 
satisfies 



\^-ip h \2<C x inf \ 1 p-X h \ 2 +C 2 y/\^h\\^-^ H \i, 
\ h ex h 



(6.4) 



where d = 2 + Re Rq- 1 L 2 + Re 2 Ro- 1 T 1 \\F\\_ 2 and C 2 = 2Re 2 Ro~ 1 T l C ||F||_ 2 . 



Proof. Subtracting (5.2) from (3.1) and letting \ — X h G C X yields the error 
equation: 

a(V - ^, x 71 ) + 4>> x h ) - b(ip H ; ip\x h ) + <i> - Tp h ,x h ) = o, v x ' 1 e X h . (6.5) 

Now, adding the terms 

-b(iP;ip h , X h ) + b(ip;ip h ,X h ) 



to (6.5) gives 

- V", x' 1 ) + 6(V; V> - i> h , x h ) + b(i/> -ip H ;i> h , x h ) + c{4> - ip h ,x h ) = o, V x h € x h . 

(6.6) 

Take \ h e X h arbitrary and define e := ip — ip h = rj — $ h , where <& h = ip h — \ h and 
i] = ip — A . Equation (6.6) becomes 

a($ h , X h ) + b(i>;$ h ,X h ) + c($ h ,x h ) = 

a(r,, X h ) + b(ij;ri,X h ) + b(ip-ip H ;i> h ,x h ) + c(r 1 ,X h ), V/eI k . 

Since ( |6.7[ ) holds for any x h G Jf' 1 it holds in particular for \ h — ® h G -^"'\ which 
implies 

a($'\ $ h ) + b(i/r, $\ $' 1 ) + c($'\ = 

a{n, <f> h ) + 6(^; r), <S> h ) + b{iP - ip H ; ^ h , $ h ) + c( V , $ h ). 



Note that c(ip,x) — — c(x>V0 which implies c(3> h ,3> h ) — 0. Also, b(ip;x>x) = 
and so (6.8) becomes 

a($ h , $ h ) = a(?7, + b(ip; f], <f> h ) + 6(V - ip h ; ip h , $ h ) + c{q, $ h ). (6.9) 



Lemma 7 



yields 



Now rewriting b(ip — tp h ; ip h , using 

a($\ $ h ) = afa, $ /l ) + b(ip; n, $'*) 

+ 6 (^ h ; ® h , i> - i> h ) + b ($ h ;i> h , i> H - if}) + c(n, <S> h ). 
Using the error bounds given in Lemmas [l] [2] |4j and [6] in (6.10) gives 
Re-^H < Re- 1 ^ \§ h \ 2 + V, |^| a |^| 2 |$ h | 2 

+ 2Ti C l^b |$"| 2 IV " >/| ln(/i) | +JJO- 1 r 2 |77| 2 
= (Eo" 1 L 2 + i^e- 1 + Lx |V| 2 ) h| 2 |$ fc | a 
+ 2Ti C |^| 2 |$ fc | 2 |^ - ^liv^M- 
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(6.10) 



Simplifying by \4> h \ 2 and using the stability estimates (3.7) in Lemma 2 and (4.2 1 in 
|Lemma"3| gives 



|$ fc |a < (l + ReRo- 1 T 2 + Re 2 Ro- 1 T 1 \\F\\^ 2 ) \r]\ 2 
+ 2Re 2 Ro- 1 T 1 C \\F\\- 2 |V - *l> H \iy/\ MW- 



(6.11) 



Adding 1 77] 2 to both sides of ( |6. 1 1 ) and using the triangle inequality gives 

< (2 + ReRo- 1 r 2 + Re 2 Ro- 1 r 1 \\F\\_ 2 ) |r?| 2 
+ 2i?e 2 i?o- 1 r 1 C ||F||_ 2 |V - V ff |iV|ln(/i)|. 



Taking the infimum over A h e X 71 in (|6.12[) yields the estimate flO|). 



(6.12) 



□ 



In what follows, we consider both X' 1 and X ff Argyris FE spaces. We emphasize, 
however, that both | Algorithm i] and the error estimate in |Theorcm 2| remain valid for 
other conforming FE spaces, e.g. the Bell element, the Hsieh-Clough- Tocher element, 
or the Bogner-Fox-Schmit element. 

For the Argyris triangle; we have the following inequalities, which follow from 
approximation theory [2] and Theorem 6.1.1 in [8]: 



IV - < Ch 6 - j , 

\iP-tP h \j<CH 6 -i, 



(6.13) 



where j = 0, 1, 2 and ip, the solution of (3.1 1, is assumed to satisfy ip G H 6 (£1)<~)Hq (fl). 
Corollary 1. Let X h . X H C Hq (SI) be Argyris finite elements. Thenip h , the solution 
of the two-level algorithm ^Algorithm 1 ) satisfies the following error estimate: 



IV - ^% < c,h 4 + c 2 ^/\Hh)\H 5 



Proof. This follows directly by substituting the inequalities (6.13) into (6.4) 



(6.14) 



□ 



7. Numerical Results. The goal of this section is two-fold: first, we illustrate 
the computational efficiency of the two-level method, and second, we verify the the- 



oretical rates of convergence developed in Section 6 To illustrate the computational 
efficiency of the two-level method, we compare solution times for the full nonlinear 
one-level problem and for the two-level method applied to the SQGE. We choose 
coarse mesh/fine mesh pairs such that the ratio is 1/2. To verify the theoretical rates 
of convergence, we compare the theoretical error estimates to the observed rates of 
convergence from our numerical tests. For the one-level problem we rely on our orig- 
inal code that was benchmarked in [TI] . 



To this end we apply the two-level method to the SQGE (2.1 ) with Re = Ro = 1 
and exact solution 



ip(x, y) = (sin Airx ■ sin 2nyY 



Additionally, the homogeneous boundary conditions are ip 



(7.1) 



and the forcing 

function F corresponds to the exact solution \7.1\. These boundary conditions and 



,9V 
dn 



exact solution will be used in all the two-level tests that follow. 
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7.1. Practical Consideration. A key part of two-level algorithms is accessing 
a previous coarse mesh solution, i.e. finding the parent element given a child element. 
This step can negate any performance benefits if not implemented wisely. Indeed, 
let n be the number of elements in the FE discretization. For the unit square, a 
naive search across every element takes 0(n/2) operations. This procedure may be 
improved with a binary search, which is summarized in |Algorithm 2| 

We note that every element on the fine mesh corresponds to exactly one element 
on the coarse mesh. However, a coarse mesh element may correspond to multiple 
elements on the fine mesh. 

Algorithm 2 Given an element on the fine mesh determine the parent element on 

the coarse mesh. 

Before examining the fine mesh, sort the coarse mesh elements by their centroid values. 

Step 1: Select an element on the fine mesh and compute its centroid. 

Step 2: Perform a binary search across the coarse mesh elements until the difference 

between the x-values of the fine mesh centroid and coarse mesh centroids is 

less than H, the coarse mesh step size. There should be many elements that 

fit this condition; save them as a list. 
Step 3: Search through this list until we find the correct coarse mesh element (that 

is, the centroid of the fine-mesh element is an interior point of the correct 

coarse mesh element). 



For the considered unit square, the binary search will examine on average log(n) 
elements, while the linear search component involves at most y/n/2 elements. There- 
fore the search requires a 0(y / n/2) number of element checks. Profiling results in- 



dicate that using Algorithm 2| to identify parent elements takes much less time than 
either setting up or solving the systems, so this approach is fast enough that lookup 
time does not contribute significantly to overall solution time. 

7.2. Computational Efficiency. To illustrate the computational efficiency of 
the two-level method, we compare the simulation time for the standard one-level 
method (i.e. the full nonlinear system, without the two-level method) with the simu- 
lation time for the two-level method. 

In Table 7.1 the L 2 -norm of the error (e^), the iJ 2 -norm of the error (e#2) 
and the simulation times are listed for various mesh sizes. For each fine mesh, we 
choose a coarse mesh that ensures the same order of magnitude for the errors in the 
one-level and two-level methods. For small values of the fine mesh size, h, the two- 
level method was significantly faster than the one-level method. The errors in the 
H 2 -norm were nearly identical, while the error in the L 2 -norm were generally of the 
same order of magnitude. We also note that the tolerance in Newton's method seems 
to cause a plateau in the L 2 -norm of the error. The results in Table 7.1 are illustrated 



graphically in Figure 7.1 In this figure the simulation times of the one-level method 
(green) and of the two-level method (blue) are displayed for all the pairs (h, H) in 
Table 7.1 Figure 7.1 clearly shows that as the number of degrees of freedom (DoFs) 



increases, the computational efficiency of the two-level method increases as well. 

7.3. Rates of Convergence. The goal of this subsection is to verify, numer- 
ically, the theoretical rates of convergence in | |6.14| | of Corollary]!] Unlike the the- 
oretical error estimates for the one-level method developed in 14J, for the two- level 
method we must verify rates of convergence for two different meshes: the fine mesh, 
h, and the coarse mesh, H . 



H 


h 


DoFs, H 


DoFs, h 


e L 2 


e H 2 


time, s 


— 


0.05146 


— 


4362 


4.286 x 10~ 8 


1.648 x 10 


-3 


3.328 


0.1083 


0.05146 
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1.092 x 10 _r 


1.709 x 10 


-3 


2.372 
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0.02561 
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16926 


5.748 x 10~ 10 


1.009 x 10 


-4 
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0.02561 
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1.016 x 10 


-4 
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4.751 x 10 -11 
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-5 
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-6 


59.03 
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0.009659 
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2.382 x 10 
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0.02035 
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95.93 
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0.006854 
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0.01436 


0.006854 
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0.006374 
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3.912 x 10 -11 


3.846 x 10 
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0.01277 


0.006374 
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2.309 x 10 -11 
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0.005264 
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-7 
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0.01101 


0.005264 
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389994 


6.495 x 10~ n 


2.087 x 10 


-7 


397.7 



of 



Table 7.1 

Comparison of one-level and two-level methods: the L 2 -norm of the error (e L 2 ), the H 2 -norm 
the error (e H 2 ) and simulation times. 



- One-level Algorithm 
" Two-level Algorithm 




1 1.5 2 2.5 

Degrees of Freedom (DoFs) 



Fig. 7.1. The simulation times of the one-level method (green) and of the two-level method 
(blue) are displayed for all the pairs (h,H) in \Table 7.i\ 



To verify, numerically, the theoretical rate of convergence with respect to H, given 



in (6.14), we fix h to a small value and we vary H. Thus, the total error in ( 6.14| will 
be dominated by the H term, i.e. the total error will of order 0(H 5 ). In Table 7.2 
we fix h — 0.0063 and we vary H. The error in the L 2 -norm (e^), the error in the 
i/ 2 -norm (e#2), and the rate of convergence with respect to H are listed in Table 7.2 
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h 


DoFs, H 


DoFs, h 


e L 2 


e H 2 


H 2 order 
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0.0063 


350 
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3 


4.54 x 10" 


i 
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0.0063 
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4838 


296710 


1.86 x 10" 


-7 


4.92 x 10" 


4 


5.04 


0.05 


0.0063 


18886 


296710 


2.11 x 10" 


-9 


1.32 x 10- 


-5 


5.2 


0.025 


0.0063 


74630 


296710 


2.. 18 x 10" 


-9 


6.02 x 10" 


-7 


4.45 



Table 7.2 

Two-level method: the L 2 -norm of the error (e L i ), the H 2 -norm of the error (e H 2 ), and the 
convergence rate with respect to H 



The rate of convergence follows the theoretical rate predicted in (6.14) (i.e. fifth- 
order). For the last mesh pair, however, the rate of convergence appears to drop off. 



This occurs because, for small values of H, the total error in (6.14) is not dominated 
anymore by the H term. 

To verify, numerically, the theoretical rate of convergence with respect to h, given 
in (6.14), we must proceed with caution. The reason is that a straightforward ap- 
proach would fix H and let h go to zero. This approach, however, would fail, since 
the H term would eventually dominate the total error. To avoid this, we consider the 
following scaling between the mesh sizes: 



H = C h, 



(7.2) 



where C > 1. The scaling in (7.2) implies that the total error in (6.14) is of order 
0(h 4 ). Indeed, the second term on the right hand side of (6.14) now scales as follows: 



C 2 y/\Hh)\H 5 « C 2 CV\Hh)\h 5 
« 0(ft 4 ), 



(7.3) 



where in the last relation in ( |7.3| ) we used the fact that yf\ m(/i)| h — » when h — » 
(which follows from 1'HospitaPs rule). 

Remark 1. We emphasize that the scaling in (7.2) is not needed in the two-level 
algorithm. We only use it in this subsection to monitor the convergence rate with 
respect to h. 

In this subsection, we consider C = 2 in (7.2). We note, however, that any other 



constant C > 1 could be used in (7.2). With this choice, we are now ready to verify, 
numerically, the theoretical rate of convergence with respect to h given in (6.14), 

In Table 7.3 for various mesh size 



which, as shown in ( |7.3[ ), will be of order 0(h 4 
pairs (H = 2h,h), we list the L 2 -norm of the error (e^), the i? 2 -norm of the error 
{&h 2 )i an d the rate of convergence. The rate of convergence follows the theoretical 
rate predicted in (6.14) (i.e. fourth-order). 



8. Conclusions. In this paper, we proposed a two-level FE discretization of the 
(nonlinear) stationary quasi-geostrophic equations. The two-level algorithm consists 
of two steps. In the first step, the nonlinear system is solved on a coarse mesh. In 
the second step, the nonlinear system is linearized around the approximation found 
in the first step, and the resulting linear system is solved on the fine mesh. 

Rigorous error estimates for the two-level FE discretization were derived. These 
estimates are optimal in the following sense: for an appropriately chosen scaling 
between the coarse mesh, H , and the fine mesh, h, the error in the two-level method 
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H 


h 


DoFs, H 


DoFs, h 


e L 2 


e H 2 


H' 2 order 


1.1 


0.43 


38 


106 


1.31 x 10" 


-2 


5.832 x 10 u 




0.43 


0.21 


106 


350 


2.25 x 10" 


-3 


6.61 x 10" 1 


2.56 


0.21 


0.10 


350 


1270 


4.40 x 10" 


-5 


3.65 x 10~ 2 


3.96 


0.10 


0.05 


1270 


4838 


1.90 x 10" 


-7 


2.00 x 10~ 3 


4.10 


0.050 


0.025 


4838 


18886 


8.95 x 10- 


10 


1.20 x 10~ 4 


4.04 


0.025 


0.013 


18886 


74630 


1.36 x 10- 


10 


7.40 x 10~ 6 


4.02 


0.013 


0.0063 


74630 


296710 


2.22 x 10" 


-9 


4.68 x 10~ 7 


3.99 



Table 7.3 

Two-level method: the L 2 -norm of the error (e L 2 ), the H 2 -norm of the error (e H 2 ), and the 
convergence rate with respect to h 



is of the same order as the error in the standard one-level method (i.e. solving the 
nonlinear system directly on the fine mesh). 

Numerical experiments for the two-level algorithm, with the Argyris element, 
were also carried out. The numerical results verified, numerically, the theoretical 
error estimates, both with respect to the coarse mesh size, H, and the fine mesh 
size, h. Furthermore, the numerical results showed that, for an appropriate scaling 
between the coarse and fine meshes, the two-level method significantly decreases the 
computational time of the standard one-level method. 
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